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Summary. — 

Using the Gaussian Ansatz for the monomer-monomer correlation functions we de- 
rive a set of the self-consistent equations for determination of the conformational 
state in the bead-and-spring copolymer model. The latter is based on the Edwards 
type effective free energy functional with arbitrary two-body interaction matrix. 
The rate of conformational changes in kinetics may be expressed via the instanta- 
neous gradients of the variational free energy functional in the space of the averaged 
dynamical variables. We study the equilibrium and kinetics for some periodic and 
random aperiodic amphiphilic sequences in infinitely diluted solution. Typical equi- 
librium phase diagrams are elucidated and the conformational structure of different 
states is discussed. The kinetics of compaction of an amphiphilic copolymer to 
the globular state proceeds through formation of locally frustrated non-equilibrium 
structures. This leads to a rather complicated multistep kinetic process. We ob- 
serve that even a small modification in the primary sequence of a copolymer may 
significantly change its kinetic folding properties. 

PACS 36.20.-r - Macromolecules and polymer molecules. 
PACS 64.60.My - Metastable phases. 
PACS 64.70.Pf - Glass transitions. 
PACS 01.30.Cc - Conference proceedings. 



1. Introduction 



Conformational transitions of heteropolymers in dilute solutions attract special at- 
tention of theorists in recent years [1]. There is a significant body of works carried out on 
concentrated solutions, mixtures and blends of copolymers based on the density variables 
formulation [2]. The infinite dilution limit, for which such approaches cannot be used, 
is nevertheless very important for understanding the fundamental interactions of macro- 
molecules not beset with the complications due to the aggregation phenomena. Apart 
from a purely academic interest, the single chain problem appears to be of paramount 
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importance for building up increasingly sophisticated models of biopolymers — the goal 
that seemed too remote only a decade ago. The progress in the modern biotechnology 
is often impeded by the inability to solve a few fundamental theoretical problems among 
which there is one grand challenge — given the primary sequence of macromolecule 
find its conformational structure at a given equilibrium, steady, or more interestingly, 
noncquilibrium state. The next part of the puzzle would be to find the relation between 
the conformational structure and biofunctionality of say a protein, but this sometimes 
can be explained by location of certain active sites on the surface of a protein globule. 

The computational difficulties are clearly immense for even not so long protein 
chain of just about 60 amino acid residues. Therefore at the moment one hopes to 
develop merely some oversimplified models of heteropolymers abstracting from the ac- 
curate molecular level description. Such models are obtained by coarse-graining in a 
fashion normal for statistical mechanics by integrating out many of the local degrees 
of freedom. We hope to be able to deduce the mesostructure of a macromolecule and 
to obtain many of the important observables. The approach that we employ here is 
in its essence a nonequilibrium extension of the Gibbs-Bogoliubov variational principle. 
The main strength of the Gaussian self-consistent (GSC) method, which we [3, 4, 5] 
and others [6] were developing first for the simple homopolymer and in Ref. [7] we have 
brought to the most general form capable of describing copolymers, is in that it produces 
the complete set of mean squared distances between pairs of monomers and thus the 
conformational structure. 

In this work we do not concentrate on the formal derivation of the GSC equations. 
This is presented in some detail in Ref. [7]. First, we shall introduce a crucially important 
modification to the model itself by adding a new so-called self-interaction term. We prove 
that without this regularising term the theory is plagued with divergences if at least one 
second virial coefficient is negative. We show that although the new term is seemingly 
negligible for long chain lengths, and indeed is not required for the coil, it does in fact 
the trick of correcting the structure of the dense globular state. In the rest of the work 
we carry out extensive study of the equilibrium properties and the folding kinetics for 
a number of examples of hctcropolymer sequences. We also discuss how the spin glass 
behaviour arises for sufficiently random copolymer sequences. 

2. The Model and the GSC Equations 

The main variables in the coarse-grained description of the polymer chain [8, 9, 10] 
are the spatial monomer coordinates X„, where n is the monomer number. The solvent 
molecules are excluded from the consideration by integrating out their degrees of freedom 
from the path integral representation for the partition function. The resulting monomer 
interactions are represented by the effective free energy functional (EFEF), 



where for heteropolymers u^Jy are in principle allowed to have any dependence on the 
site indices {n} = {n\, . . . , nj}. The first term in Eq. (1) describes the connectivity of 
the chain with I called the statistical segment length. There are also volume interactions 
represented by the virial-type expansion [9, 10] in Eq. (1). They reflect the hard-core 
repulsion, the weak attraction between monomers and the effective interaction mediated 
by the solvent-monomer couplings. 
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Here we consider the following choice of site-dependent second virial coefficients in 
Eq. (1), 

(2) tt W =fi ( a)+A fn+£2l | 



and ui J \ — u (J ^ for J > 2. This corresponds to the case of amphiphilic heteropolymers, 
for which monomers differ only in the monomer-solvent coupling constants. The mean 
second virial coefficient, rS 2 \ is associated with the quality of the solvent and the pa- 
rameter A is called the degree of amphiphilicity of the chain. The set {cr n } expresses 
the chemical composition, or the primary sequence of a heteropolymer. Here we consider 
the case when variables a n can take only three values: —1,1 and corresponding to 
hydrophobic V, hydrophilic '6' and "neutral" 'c' monomers respectively. 

It is assumed that the long timescale evolution of the conformational state is well 
represented by the Langevin equation, which upon neglecting the hydrodynamics may 
be written as, 

* d ^ 9H 

where (b is the "bare" friction constant per monomer and the Gaussian noise, rj„, is 
characterised by the second momentum, 

(4) {^{t)^ l {t'))^2k B TC b 8 a - a '5 n ^5{t-t'), 

where the Greek indices denote the spatial components of 3-d vectors. 

The main idea of the Gaussian self-consistent method is to choose the trial Hamil- 
tonian, H 0} as a most generic quadratic form, with matrix coefficients depending on 
time, 

(5) H {t) = ^V nnl {t)X n {t)X nl {t). 

nn' 

This corresponds to the Gaussian distribution of the inter-monomer distances, (X m — 
X m /) 2 . Thus, the two-body monomer-monomer correlation function, that is the prob- 
ability density to find the monomer m' at the distance r from the monomer m, will be 
given by, 

(6) h%> m ,(r; t) ee (S(v X m + X ro ,)> = /0 _ n 1 y^m ex P 



(27r£> rom ,(i))3/2 2D w (t) 
Here D mm > is the matrix of the mean squared distances between monomers, 

(7) D mm ,{t) ee i ((X m (t) - X m ,(t)) 2 ) . 

Obviously, choosing Eq. (5) as the trial Hamiltonian is equivalent to replacing the non- 
linear Langevin equation (3) by a linear stochastic ensemble, 

(8) tbjXn = -£ V nn ,(t) X n , + Vn (t). 



4 



YU.A. KUZNETSOV , E.G. TIMOSHENKO 



The time-dependent coefficients are chosen at each moment in time according to the 
criterion, 

where (. . .} denotes the averaging over the trial ensemble. At equilibrium these equa- 
tions become exactly the extrema conditions for the trial free energy in the Gibbs- 
Bogoliubov variational principle based on minimising the variational free energy, A = 
-k B T\ogTrcxp(-H /k B T) + (H - H ) Q , with respect to V nn >. 

For details of calculations we refer the reader to Refs. [5, 7]. Here we present the 
final form of the kinetic GSC equations which describe the time evolution of the mean 
squared distances between monomers (7). It turns out that the equations can be written 
in terms of instantaneous gradients of the variational free energy, A = £ — TS, 

(10) %iD mm ,(t) = -lY(D mm ,,(t)-D m , m ,,(t))f (>A <>A 



2 dt 3^ \dD mm n{t) dD m > m n{t) 

m" x 

The energetic and the entropic contributions in the free energy can be completely ex- 
pressed in terms of the mean squared distances D mm > (t) , 

oo (J) 

(ID *=^|nE + EE fi^gW detA^))- 3 / 2 + 

n J=2 {n}> V ' 

(12) S = ^fcsfogdcti?^- 1 ), R nn , = ^2 E D nm,n'm', 

mm' 

where we have introduced the four-point correlation function and the matrix A^ -1 ), 

(13) D mrn ' nn ' — —{Dm'n ~t~ D mn ' Dmn ^m'n')' 

(14) A^ — D nini+lt „ inj+1 . 

In Eq. (12) we have the determinant of the truncated matrix to exclude the zero 

eigenvalue related to the translational invariance for the centre-of-mass of the system. 
In the second term in Eq. (11), which is responsible for the volume interactions, the 
summation is taken over not coinciding indices, ni / n 2 ^ . . . ^ nj. 

Before proceeding with further discussions of the GSC equations let us introduce 
some observables. These include the mean squared radius of gyration, 

(15) ^ = ^E^'' 

mm' 

and the micro-phase separation (MPS) order parameter, 

(16) * = ^E^ A ^ 1 ^ (^) 2 = ^E^- 

9 mm' m 

The MPS parameter describes the degree of correlation between matrices of the relative 
two-body interaction, (<r m + <r m /)/2, and the mean squared distances, D mm >. 
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3. The Self interaction Energy Term 

The appearance of the last term in Eq. (11), £ S i, is somewhat more nontrivial. In 
fact, in the EFEF of the model (1) we have discarded terms with two or more coinciding 
indices in the three- and higher body contributions. These terms come formally from the 
virial-type expansion, but each of them gives a singular contribution to the mean energy 
(11). It turns out, however, that upon suppressing these terms there appear additional 
pathological solutions of the GSC equations with singular free energy if at least one 

(2) 

element of the two-body interaction matrix, u ' ,, becomes negative. This is easy to 
see. Indeed, consider volume interactions of just three monomers under condition that 
the mean squared distances from monomers '0' and '1' to '2' are equal to each other, 
-Do, 2 = -Di,2 = D. These interactions produce the mean energy contribution, 

h7 , c + , 1 ( (2) , 6u^(2n)-^ \ 

[ ' (27^)3/2 + (2tt j Do,i) 3/2 V 0,1 (D- A>,i/4) 3 / 2 J 

(2) 

In the case when Uq { < and monomer '2' is placed away from monomers '0' and '1', D > 

A),i/4+ (6m( 3 )(2tt)- 3 / 2 /|m^|) 2 / 3 , obviously in the limit £> ,i -» the energy possesses 
a singular minimum, £3 — > —00. As for the free energy the logarithmic divergence of the 
entropy could not change the situation, thus A — > —00 as well. One can show that the 
inclusion of more monomers in the chain or of higher than three-body interactions does 
not improve the situation, but produces more and more of such pathological solutions. 
The reason we have not discussed this problem in our previous considerations is that we 
have accounted for the additional symmetry properties of monomer-monomer distances, 
which come from the symmetry of the EFEF (1). For example, in the case of the ring 
homopolymer, due to the inverse symmetry [5], we assumed that for any indices m, 
to' the following mean squared distances are equal, D m ^ m > = D m ,2m-m' = D2m>-m,m> ■ 
This provides sufficient repulsion coming from three-body term to preclude pathological 
solutions. 

Thus, in a more general case, where no symmetry properties could be assumed for an 
arbitrary sequence, the standard procedure of suppressing terms with coinciding indices 
is not satisfactory. Fortunately, it could be remedied by using another prescription — 
replacing the terms with coinciding indices by the so-called self-interaction terms. Here 
we propose the prescription for three-body interaction which is sufficient for our current 
purposes, 

(18) 4 3) = c 3 u^ ]T /<5(X ro X ro ,)\ = c 3 u^ ]T D~l„ 

where C3 = 3 is a combinatorial factor related to the three possible ways of having 
coinciding pairs of indices in a triple summation. Obviously, the higher negative power 
of D mm i in (18) compared to the two-body term in (11) prevents one monomer from 
falling on another. 

It would be interesting to consider the coefficient C3 in Eq. (18) as an independent 
parameter and to discuss how the inclusion of this term would change the equilibrium 
and kinetics for the ring homopolymer [11]. Note that for the ring homopolymer we can 
reduce the number of independent elements in the matrix of mean squared distances, 
D m m', by the factor of N, since due to the translational symmetry, D mtTn > = D m+ k, m '+k 
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for any k and m, m' . This symmetry allows us to reduce significantly the computational 
requirements. 

In the good solvent regime, > > 0, we find no substantial change caused by 
the self-interaction term (18). The deviations in the mean squared distances and in the 
squared radius of gyration typically are less than 1% even for short chains. This is only 
natural. In this regime this term is subdominant for long chains and can be neglected at 
all in the thermodynamic limit. 

However, the size and the structure of the homopolymer globule, which exists in the 
region < 0, can change significantly. In Fig. 1 we exhibit the diagram of globular 
states. One can see that actually two different globular states are possible, which we 
call the non-compact globule (as one can see from Table I this globular state possesses 
a higher value of R 2 g ) and the liquid-like globule. The former is the thermodynamically 
stable state at comparatively small c 3 , whilst the latter exists at higher values of C3 and in 
the region of large . For large negative the transition between two globular states 
becomes discontinuous and the transition line goes nearly vertically at approximately 
c 3 « 1/2. However, the collapse transition from the extended polymer coil to the liquid- 
like globule remains second order. We found that this phase diagram is quite independent 
of the degree of polymerization, at least in the range we could study numerically, 30 < 
N < 200. Thus, introducing of the term (18), which is subdominant for large values of 
N compared to other terms in (11), nevertheless dramatically changes the properties of 
the globular state even for very long chains. 

Now let us compare the structure of the non-compact and the liquid-like globules. In 
Fig. 2 in the left-hand side we present the mean squared distances, which are symmetrical 
with respect to the line m = N/2, i.e. £>o,|m-ro'| = D o,N-\m-m'\i f° r values of C3 = 0, 3. 
In the right-hand side of the figure we draw the same quantity obtained from Monte 
Carlo simulations for a ring chain of hard spheres with the Lennard- Jones attraction 
[15]. Since the parameters of excluded volume interactions in the model used for Monte 
Carlo simulations are expressed in different terms, only shape and scaling behaviour, but 
not the absolute values of each D mm >, should be compared. We see a quite remarkable 
agreement here. This is a strong argument in favour of the liquid-like globule. Its mean 
squared distances have a typical saturation regime after comparatively small m. This is 
known to be true of the globule from Monte Carlo simulations and previous results from 
other methods [12]. The function, Do m for the non-compact globule has much more 
convex shape which, in fact, reflects the effective monomer repulsion on large distances 
(see Eq. (17). Also, from Table I we can see that the asymptotic scaling in the degree 
of polymerisation for the size of the liquid-like globule, Rg ~ iV 2 / 3 , is reached starting 
from sufficiently small N in agreement with lattice Monte Carlo simulations [13], whilst 
from much larger N for the non-compact globule [4] . 

Now, let us turn our attention to the kinetics at the collapse transition, after an 
instantaneous change of the value of the second virial coefficient from a positive to a 
negative value, the latter corresponding to the globular equilibrium state. In Fig. 3 we 
exhibit the evolution of the mean squared radii of gyration after quenches to the liquid- 
like and the non-compact globule. We found that the early stage proceeds in a similar 
manner, while the middle or "coarsening" kinetic stage is somewhat slower for quenches 
to the liquid-like globule, but the final stages here are much faster. This can be seen 
from Tab. I, where we present the values of the total collapse time [14] and the final 
relaxation time, 17 for different sizes, N. The most striking thing here is that the effective 
exponent of the final relaxation time is much smaller than what we have earlier expected, 
7/ = 5/3 (see Ref. [3, 4]). This also shifts down the effective exponent of the "total" 
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collapse time, since the latter is a cross-over between the exponent of the coarsening 
stage, -f m = 2, which is unaffected by the self-interaction term, and the exponent of 
the final relaxation, jf. In principle, the final relaxation time and its exponent in the 
degree of polymerization can be determined without appealing to numerical solution of 
the kinetic GSC equations. From Ref. [3] we have, r/ ~ iV.Fi, where Fi is the first 
normal mode for the final equilibrium globular state, 



J_^ ( ,,.^l/) :: 



A good approximation for the function D^ m in the liquid-like globule would be the 
following Dom is linearly increasing function of m until some value Do mi = D, where the 
TV-dependence is, D ~ m\ ~ TV 2 / 3 and then it remains constant. This approximation 
can also be obtained from the Lifshitz theory [12]. Now one can see that T\ ~ iV 1 / 3 and, 
correspondingly, jf = 4/3, which is close to the value of the effective exponent in Tab. I. 
Thus, on the contrary to the coil state, for the liquid-like globule the first normal mode, 
T\ , neither gives the main contribution to the mean squared radius of gyration, nor even 
has the same scaling law in the degree of polymerization, N. 

Inclusion of the self-interaction term (18) also strongly affects the phase diagram and 
kinetic properties of rigid chains. We shall consider these questions in a separate work 
[15]. Such modifications arc rather welcome and allows one to make the GSC method 
more accurate for the dense globular states, where its validity has been less established. 



4. Equilibrium Properties of Copolymer Sequences 

The GSC equations (10) have been studied numerically using the fifth order Runge- 
Kutta algorithm with adaptive time step [16]. In fact, for copolymer sequences, consid- 
ered in the current work, the time step during numerical integrations of Eqs. (10) varies 
approximately 100 times. Thus, any integration scheme with fixed time step is rather 
unreliable here. 

In studying the equilibrium we consider only stationary points of Eqs. (10), i.e. the 
limit t — > oo. If for some set of interaction parameters, and A, one obtains several 
stationary states, one should compare the values of the variational free energy, A. The 
deepest minimum of the free energy corresponds to the thermodynamically stable state, 
the rest of the solutions to metastable states. Here we consider copolymers with the ring 
topology, though the current treatment may be easily extended for study of copolymers 
with any other topology just by changing the spring term in Eqs. (1, 11). 

Typical phase diagrams in terms of the mean second virial coefficient, vS 2 \ and the 
amphiphilicity, A, in Eq. (2) for some "random" and periodic sequences arc presented 
in Figs. 4-6. In the region > and for small values of amphiphilicity, A < 5, typical 
conformations of copolymers are akin to the homopolymer extended coil. By decreasing 
it*- 2 ) to the negative region the chain undergoes the continuous collapse transition, simi- 
larly to what we observed in Sec. 3. The collapse transition is characterised by a rapid 
fall of the radius of gyration, i? 2 , (15) and the change of the fractal dimension, v (see 
Tab. I). 

The collapse transition for larger values of amphiphilicity turns out to be more com- 
plicated, and essentially dependent on the sequence. The globular state for large values of 
A is different from the liquid-like globule. It is characterized by somewhat higher value 
of the radius of gyration and extremely large value of the MPS order parameter (16), thus 
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we call this state the MPS globule. The MPS globule is separated from the liquid-like 
one by a weak continuous transition (see Figs. 4-6). In the case of long blocks (Fig. 6) 
the collapse transition to the MPS globule becomes discontinuous (first- order-like) . The 
spinodals F and I" designate the region where two distinct states corresponding to the 
coil and the MPS globule can be found. The depths of the free energy minima become 
exactly equal on the transition curve I in Fig. 6. 

However, for a wide class of sequences, for example for aperiodic sequences in Figs. 
4, 5, the phase diagram at large amphiphilicity, A, is much more complicated. Starting 
from some value of A in some intermediate region of there appear additional solutions 
corresponding to local minima of the free energy. The broad region where this could take 
place is bounded by the curves F and II" in Figs. 4, 5. With increasing A the number of 
such solutions grows quickly. Significantly, in the region of the phase diagram, between 
curves I and II in Figs. 4, 5, some of these possess the lowest free energy value, thus 
being the thermodynamically stable state. Since the number of such solutions is rather 
high even for short sequences and their number grows quickly with the chain length, 
we do not attempt to draw all their boundaries of (meta)stability. We shall call them 
collectively as frustrated phases, explaining this terminology below. 

Now let us compare the phase diagrams in Figs. 4 and 5, the latter corresponding 
to the sequence twice longer than the for former. An interesting observation is that 
the region between spinodals F and II", designating where the frustrated phases can 
exist, expands dramatically with increasing chain length. The same concerns the region 
of thermodynamically stable frustrated phases between curves I and II. More exactly, 
these regions expand downwards and to the left, so that the position of curves I' and I 
change slightly with increasing N , whilst curves II and II" depend significantly on the 
size of the system. For rather long chains we may expect that the regions of stability 
and metastability of the frustrated phases will continue to expand downwards and to 
the left, so that the lines II and II" will become nearly vertical, displacing the region of 
stable MPS globule. Probably, for most of long heteropolymer chains the MPS globule 
does not exist as thermodynamically stable state, becoming stable only for some special 
sequences. Unfortunately, we can not proceed with numerical solution for much larger 
system sizes, N, since the calculational expenses grow in TV as N 3 per iteration and also 
the total number of frustrated states becomes huge for large system sizes. This diversity 
and a special foliating structure of various branches leads in the thermodynamic limit to 
what is known as a spin glass frozen phase [17] of random copolymers. 

Let us consider the conformational structure of the frustrated states for the copoly- 
mer consisting of repeating 'a6' blocks. The phase diagram of this sequence also exhibits 
the thermodynamically stable frustrated phases [7] starting from approximately N = 28. 
In Figs. 7 and 8 we present the dependence of the observables in (15, 16) at a quasistatic 
change of the mean second virial coefficient, tS 2 \ from the coil state to the MPS globule 
and back. From these pictures one can see in the intermediate region only a few (seven 
to be precise) of all possible frustrated phases. The values of the radius of gyration and 
the MPS order parameter are intermediate for these solutions, lying between those of 
the coil and the MPS globule. In this sense, we can call them non-fully compacted and 
misfolded states. 

In the series of pictures in Fig. 9 we present the matrix of mean squared distances, 
D m m>, for the copolymer, consisting of "aV blocks at some values of the mean second 
virial coefficient, u^ 2 \ The set of D mm > in the GSC method completely determines the 
conformational structure of any, equilibrium or kinetic, state. For positive (see 
Fig. 9a) the mean squared distances possess the structure typical for the extended coil. 
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The elements of the matrix, D mrn >, increase monotonically on moving away from its 
diagonals towards the distance of half-ring along the chain. Thus, the D mm i matrix 
may be approximated here by a monotonically increasing function of the chain distance, 
m — m'\. Decreasing the mean second virial coefficient, u^ 2 \ causes the copolymer to 
pass through frustrated states, in Figs. 9b-9d, finally reaching the MPS globule state (see 
Fig. 9e). The characteristic feature of the D mm i matrix in a frustrated state is that it 
possesses some number of monomer groups having smaller distances between each other 
than between monomers from other groups. Clearly, such a group represents a cluster 
of monomers, so that here the copolymer chain forms a set of clusters (approximately 
8 clusters in Fig. 9b, 4 clusters in Fig. 9c and 2 clusters in Fig. 9d), each consisting of 
monomers nearest along the chain. 

The internal structure of each cluster is similar to the structure of the final MPS 
globule. For the simplest copolymer sequence presented here, the D mm i matrix away from 
the diagonal can be divided into sub-matrices of the size 2x2 (see top of the Fig. 9e), 
each of which possesses approximately the same structure: small values in the upper 
left corner correspond to the mean squared distances between hydrophobic monomers; 
large values in the lower right corner, which are the mean squared distances between 
hydrophilic monomers; and off-diagonal elements are nearly equal and correspond to 
the distances between different species. Thus, for the MPS globule, such structure of the 

(2) 

Dmm' matrix reflects the structure of the two-body interaction matrix, u mm , . Obviously, 
the higher correlation between these two matrices is manifested in the higher value of 
the MPS order parameter. For another, more complicated, sequences the D mm i matrix 
in the MPS globular state has a different structure, but still it resembles in some way 
the interaction matrix, u^ m , away from the diagonal, being distorted by the spring 
interactions in the elements close to the diagonal. Note that for the liquid-like globule 
the pattern of the D mm i matrix looks sufficiently trivial. Namely, it has narrow band of 
dark color along the main diagonal with its intensity quickly decreasing, as one can see 
from Fig. 9f. 

From Fig. 9 one can see that the interaction symmetries imposed by the EFEF (1) are 
spontaneously broken in the region of the phase diagram corresponding to the frustrated 
phases. Obviously, in our case the number of possible ways to break the symmetry 
is enormously huge and moreover it grows exponentially with increasing system size. 
Despite the kinematic symmetries are not present for arbitrary sequences, the structure of 
the phase diagram (see Figs. 4, 5) and behaviour of main observables remain very similar. 
It is the particular structure of D mm i matrix, number and the shape of boundaries of 
frustrated phases that are quite sensitive to the sequence. The symmetry that may 
be broken in this case has a subtler meaning and may be expressed in terms of the 
replica formalism [17]. Consider, for example, two nearly identical blocks with nearly 
identical surroundings. In the coil and the MPS globular states one should expect the 
conformational matrices of these blocks to be nearly equal to each other. On the contrary, 
the small difference in the interactions in the region of frustrated phases may lead to a 
huge difference in the conformations. Thus, it is by no means surprising that the replica 
symmetry breaking in the case of periodic systems takes such an explicit manifestation 
in the breaking of the block translational symmetry. 

An important point here is that the frustrated phases become dominant in an inter- 
mediate region of the phase diagram not due to a low mean energy, but mostly due to a 
higher entropy The MPS globule is entropically unfavourable there because the overall 
shrinking force is insufficiently strong. Also, in the regime of nearly compensating repul- 
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sion and attraction between monomers it is more preferable to achieve phase separation 
on a smaller, than the globular, scale by forming clusters. 

5. Folding Kinetics 

Here we shall consider the time evolution of the conformational state of a copoly- 
mer away from the initial equilibrium after it has been subjected to an instantaneous 
temperature jump that causes the two-body interaction parameters in Eq. (2), and 
A, to change. The composition of copolymer {er} and the rest of parameters remain 
the same and do not change with time. We are interested in quenches from the ho- 
mopolymer coil, where all monomers are equally hydrophilic (u^ > and A = 0), to 
the region of parameters corresponding to the MPS globular state, so that the 'a' species 
became strongly hydrophobic and the '6' species remained nearly neutral (u^ <C 
and A w |u( 2 )|). Here we consider some binary copolymer systems, consisting of the 
same number of 'a' and '6' monomers and of the same size, TV = 50, and quench, 
(u^ = 15, A = 0) -» (u^ = -35, A = 30). 

Let us first discuss the general features of the kinetics of the copolymer folding. The 
time evolution of the mean squared radius of gyration, R 2 g {t), the MPS order parameter, 
^(t), and the instantaneous free energy, A(t), is presented in Figs. 10, 11. Here lines 
A correspond to the homopolymer and serve for reference purpouses. In the case of 
the homopolymer one can see that R 2 g and A decrease monotonically to their equilibrium 
values corresponding to the liquid-like globule state, while remains identically zero. As 
for the copolymer kinetics, the first observation consists in that it proceeds much slower 
than for the homopolymer. For example, for the simplest copolymer sequence, '(60)25', 
the total collapse time [14] is more than 3 times longer than that of the homopolymer, 
other copolymer sequences collapsing even longer. The total collapse time also seems to 
be quite sensitive to the copolymer sequence. Since the micro-phase separation is one 
of the main factors in collapse of copolymers, the MPS order parameter, "J, grows in 
kinetics most of the time for all considered sequences, though rather nonmonotonically. 
Nevertheless, for some sequences it might be negative during some time in kinetics (see 
e.g. sequences E and F in Fig. lib), something we never observed at equilibrium. 

The evolution of the instantaneous free energy, A(t), depicted in Figs. 10c, 11c, is 
most unusual. Typically it proceeds through multiple accelerations and decelerations. 
The flat regions of a staircase-like function correspond to temporary kinetic arrest of the 
system in transient noncquilibrium conformations, i.e. to transient trappings of various 
members of the ensemble in their local shallow energy minima. Since such minima are 
encountered at different moments in time for different members of the ensemble, their 
influence on the overall time evolution of averaged observables is manifested in a smooth 
characteristic slowing down. Note here that the number of steps in kinetics process hardly 
can be guessed from the given primary sequence. Typically, the kinetics for sequences 
with smaller number of blocks proceeds through smaller number of steps. For example, 
folding of the periodic sequence consisting of long blocks, '(6505)5 (see line E in Fig. 11c), 
proceeds through only one, though rather long, kinctically arrested step. However, for 
the aperiodic sequence of long blocks (line F in Fig. 11c) we can see at least six such 
steps, the third one being the longest in time. 

Now let us consider the noncquilibrium conformations, given by the matrix of the 
mean squared distances between monomers, D mm >. In Fig. 12 we exhibit the D mm > ma- 
trix for '(60)25' sequence at different times in the folding kinetics. We can see that kinetic 
process proceeds through formation of locally collapsed and phase-separated clusters. 
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The initial conformation is similar to Fig. 9a, then at time t = 3.6 we can see approx- 
imately 10 clusters (Fig. 12a), at time t = 5.75 — 6 clusters (Fig. 12b) and at t = 8.6 
— 4 clusters (Fig. 12c). During later evolution these four clusters approach each other 
since, as one can see from Fig. 12d, the inter-cluster distances are much smaller than 
in Fig. 12c. Finally, the clusters unify forming the MPS globule, the D mm i matrix of 
which is quite similar to that presented in Fig. 9e. We should note here also that these 
nonequilibrium states do not possess the translational block symmetry, present for the 
system in the EFEF and the initial conditions. After some time in kinetics, t ~ 1 this 
symmetry breaks down, and restores only during final kinetic stages at time moment 
t « 15. 

However, the formal structure of the GSC kinetic equations (10) is such that they 
yield a symmetric solution at any moment in time if one proceeds from the symmetric 
initial condition. What we observe here is a spontaneous symmetry breaking in kinetics. 
Namely, at some moment in time the symmetric solution of the kinetic equations be- 
comes unstable with respect to perturbations, whether of the initial condition, or of the 

(2) 

interaction matrix, u n ^,. In the exact theory there are fluctuations that can transform 

between different spontaneously broken states in kinetics, thus destroying the unstable 

symmetric state. To describe such phenomena strictly in the framework of the Gaussian 

method, which presents an improved, but still a mean-field type of theory, we should 

include explicitly an infinitesimal symmetry breaking term e nn i to the two-body inter- 
im 

action matrix, u y m ' m ,, and consider the limit of vanishing perturbation in the solution. 
In fact, in numerical integration such a regularisation procedure is not even necessary, 
since there is always an intrinsic perturbation due to computer round off and numerical 
integration errors. Thus, if the symmetry is favourable to be kinetically broken some- 
where, numerically one obtains some spontaneously broken solutions there, rather than 
the unstable symmetric solution, unless the symmetry conditions have been imposed by 
hand. In fact, adding a rather small perturbation matrix, e nn >, changes the behaviour of 
the global observables, such as R 2 g , "J/ , A and characteristic kinetic times rather insignifi- 
cantly. However, what can be changed by including of a perturbation is the centers along 
the chain around which clusters form and grow. 

Finally, let us discuss how the folding kinetics depends on the sequence. In series of 
pictures in Fig. 10 we present kinetic processes for the simplest copolymer consisting of 
'ab' blocks (sequence B), and also for some sequences obtained by certain modifications 
in it. In the sequence denoted by C in Fig. 10 we have replaced ten short blocks by 
four of longer size. In the sequence D we have inserted only two hydrophobic and two 
hydrophilic fragments into the sequence, i.e. we have made only two permutations. For 
the former sequence the total kinetic time became approximately 1.7 times longer than 
that for the B sequence. The final values of the free energy for both sequences B and C are 
nearly the same, whilst the micro-phase separation is somewhat better for the modified 
sequence C, due to longer blocks. Thus, for that sequence the kinetic properties do not 
deteriorate very much, except that kinetic process takes essentially longer. However, the 
kinetic foldability of the D sequence is much more poor than for B and C sequences. In 
fact, here not only the kinetic process takes much longer, but the final state is different 
from the MPS globule. As one can see from Fig. 13 the final kinetic state for the D 
sequence consists of two clusters. These clusters are connected by two links, formed by 
nearly neutral fragments, '6 2 '- Further collapse of this conformation is unfavourable due 
to the entropy and partial screening of hydrophobic monomers by nearly neutral species 
'b\ Size of such misfolded state is larger and the MPS order parameter is approximately 
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twice smaller than what we could expect for the MPS globule. As we have already seen 
kinetics of sequences consisting of longer blocks is also quite sensitive to the sequence. 
Some modifications of the periodic sequence E in Fig. 11 to form aperiodic long blocks 
sequence, F, result in somewhat longer, but much more complicated, kinetic processes. 

6. Conclusion 

In this paper we applied the extended GSC method to studying the equilibrium and 
kinetic phenomena of some particular amphiphilic hctcropolymer sequences. First, we 
have revised the Gaussian theory to include the self-interaction term (18). Thus, the 
two-body term in our formalism describes the effective interactions between monomers, 
the self-interaction term is responsible for stability of the globular state, and the three- 
body interactions between distinct monomers provide the correct scaling of the size of 
the globule. The structure of the homopolymer globule is now in better agreement with 
the one obtained from other methods. We have found that the kinetic laws of the ring 
homopolymer do not change significantly due to introduction of the self-interaction term, 
except for the final relaxation time, which without hydrodynamic interaction scales with 
the degree of polymerization as, r/ ~ iV 4 / 3 . 

We have obtained typical phase diagrams for heteropolymers consisting of short and 
long blocks. For sequences with rather frustrated interactions along the chain, apart 
from the coil, the liquid-like and the micro-phase separated globular states, we have 
discovered that in a wide intermediate region of the phase diagram there may be a large 
number of frustrated partially misfolded states. Some of such states, the number of which 
grows exponentially with the chain length, become the dominant thermodynamic state 
in rather narrow domains. We may conclude that the transition to these states for large 
systems corresponds to a glassy freezing transition. 

For large amphiphilicity parameter collapse from the extended coil state to the MPS 
globule for a wide class of sequences proceeds at equilibrium through these frustrated 
phases. A typical copolymer conformation here is a set of micro-phase separated clusters. 
On approaching the MPS globule quasistatically the number of clusters decreases, reach- 
ing only one cluster in the final MPS globular state. A typical cluster size is determined 
by some characteristic range where the micro-phase separation may occur. 

We have discovered that the region of meta- and thermodynamic stability of the 
frustrated phases expands with system size displacing the region of the MPS globule. 
It is quite likely that the thermodynamically stable MPS globule for heteropolymers 
of several hundreds of monomers is possible only for a narrow class of some special 
sequences. Also it is likely that inclusion of other interactions, such as electrostatic, 
and thus modification of the two-body interaction matrix in Eq. (2), may additionally 
stabilise the MPS globule. 

We find that the kinetics of folding from the coil to the MPS globule for copolymers 
takes much longer than for the homopolymer since it is strongly affected by the presence 
of the transient frustrated states along the kinetic pathway. This leads to a complicated 
kinetic process consisting of multiple steps with pronounced slowing down and then 
acceleration in the folding rate. We also present here some preliminary studies on how 
the folding kinetics depends on the primary heteropolymer sequence. Typically, the 
kinetics for copolymers consisting of long blocks proceeds in a smaller number of steps, 
but not necessarily faster, than kinetics for short block copolymers. For the latter we 
have seen that even small modifications of the sequence may change crucially the overall 
kinetic behaviour and even the final kinetic state itself. 
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There is an interesting arena here for classification of sequences and conformational 
states. Moving in this direction would allow to develop better models for proteins, which 
sequences have been specially optimised by the biological evolution. 

Finally, we hope that the way is now open for constructing the non-Gaussian exten- 
sion to the self-consistent method. This is difficult since the closure relations for higher- 
order correlation functions are nontrivial for polymers. When that is accomplished one 
could speak of real predictive accuracy of the method. 
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Table I. - Values of the mean squared radius of gyration, R 2 g , for the globular state at equilib- 
rium, — —25 and — 10, and characteristic collapse times, rt and Tf after the quench 
u^ 1 = 15 — » —25 for different lengths of polymer ring, N . Data for cz = 3 and c 3 — corre- 
spond to the prescription with and without self-interaction term (18). The last column contains 
the effective exponent of the appropriate quantity in terms of the degree of polymerization, i. e. 
A ~ N 1 . 



N 




30 


50 


70 


100 


150 


200 


300 




exponent 


R 2 g (c 3 = 


= 0) 
= 3) 


1.19 
0.91 


1.91 
1.32 


2.59 
1.68 


3.53 
2.17 


4.96 
2.90 


6.25 
3.55 


8.52 
4.72 


2v 
2v 


= 0.86 ± 0.03 
= 0.69 ±0.01 


T t (cz = 
Tt(c 3 = 


0) 
3) 


3.9 
3.4 


10.5 
8.4 


20.3 
14.8 


40.8 
27.3 


90.7 
53.3 


160.3 
89.1 


357 
179 


It 
It 


= 1.96 ±0.01 
= 1.70 ± 0.03 


T/(C3 = 
T/(C 3 = 


0) 
3) 


1.18 
0.73 


2.97 
1.41 


5.49 
2.66 


10.2 
4.25 


20.6 
7.34 


33.0 
9.92 


65.7 
16.8 


If 
If 


= 1.74 ±0.03 
= 1.27 ±0.05 



Fig. 1. - Diagram of stability of the liquid-like versus the non-compact globules in term of the 
second virial coefficient, u^ 2 \ and the self-interaction parameter, C3. Curve I corresponds to 
discontinuous transition, curves I' and I" are spinodals. 

Fig. 2. - Plot of the mean squared distances between monomers, Dom, versus the chain index, 
m, for polymer with N = 200 for the non-compact (NCG), C3 = 0, and the liquid-like (LLG) 
globules, C3 = 3. Here — —25 and = 10. In the right-hand side we present the 
data from Monte Carlo simulations [?] for the same polymer size. Since the function Dom is 
symmetric with respect to the value N/2 we present these dependencies only on half-interval. 

Fig. 3. - Plot of the mean squared radius of gyration, R 2 g versus time, t, at the kinetics of the 
collapse transition, = 15 — » —25 to the non-compact (NCG), C3 = 0, and the liquid-like 
(LLG), C3 = 3, globules. 

Fig. 4. - The phase diagram of "random" sequence 'babca2cbac2acbzcacza2b2cac2b' in terms of 
the mean second virial coefficient, u^ 2 \ and the amphiphilicity, A. Curves (Collapse) and (MPS) 
correspond respectively to the collapse and the micro-phase separation continuous transitions. 
Curves (I) and (II) correspond to discontinuous transitions to the frustrated phases. "Spinodal" 
curves (F) and (II") bound the regions of metastability of the frustrated states. Transition 
curves and boundaries distinguishing different frustrated states are not depicted. 

Fig. 5. - The phase diagram of "random" sequence '(6o6co2c6oc2ac6 3 coc30262COC2&)2', i.e. twice 
as in Fig. 4 in terms of the mean second virial coefficient, vf- 2 \ and the amphiphilicity, A. See 
caption to Fig. 4 for more details. 

Fig. 6. - The phase diagram of "long block" copolymer sequence \bzaz)5 in terms of the mean 
second virial coefficient, v> 2 ', and the amphiphilicity, A. Since the frustrated phases here are not 
accessible, the collapse for large amphiphilicity proceeds through the discontinuous transition 
(curve (I)) and it is accompanied by micro-phase separation. Curves (F) and (I") are spinodals. 



16 



YU.A. KUZNETSOV , E.G. TIMOSHENKO 



Fig. 7. - Plot of the mean squared radius of gyration, i? 2 , versus the mean second virial 
coefficient, u' 2 ', for '(06)30' copolymer. Here and in Fig. 8 A = 30, solid lines correspond 
to values of observables in the main free energy minimum and dashed lines — in metastable 
minima. 



Fig. 8. - Plot of the parameter of micro-phase separation, <]/, versus the mean second virial 
coefficient, vS 2 \ for '(06)30' copolymer. See also caption to Fig. 7 for more details. 



Fig. 9. - Diagrams of the mean squared distances matrix, D mm i for '(06)30' copolymer and 
amphiphilicity, A = 20. Diagrams (a-e) correspond respectively to the values of the mean 
second virial coefficient, u' 2 ' = 15, —11, —21, —30 and —40. Diagram (f) corresponds to the 
homopolymer globule, for which vP^ = —25 and A = 0. Indices m, m' start counting from the 
upper left corner. Each matrix element, D mm i is denoted by a quadratic cell with varying degree 
of black colour, the darkest and the lightest cells corresponding respectively to the smallest and 
to the largest mean squared distances. The diagonal elements are not painted since by definition, 



Fig. 10. - Plot of the mean squared radius of gyration, R g (Fig. (a)), the MPS order parameter, 
$ (Fig. (b)) and the instantaneous free energy, A (Fig. (c)) versus time, t, during kinetics 
after the quench from = 15, A = to u' 2 ' = —35, A = 30 for copolymer sequences with 
N = 50. Lines A-D in the figures correspond respectively to the following sequences: 'C50' 
(homopolymer), '(60)25', '0263036203620263(60)15' and '(60)362(60)902(60)56202(60)4'. 



Fig. 11. - Plot of the mean squared radius of gyration, R 2 g (Fig. (a)), the MPS order parameter, 
* (Fig. (b)) and the instantaneous free energy, A (Fig. (c)) versus time, t, during kinetics after 
the same quench as in Fig. 10 for copolymer sequences with N = 50. Lines A, E, F in the 
figures correspond respectively to the following sequences: 'c 5 o' (homopolymer), '(6505)5' and 
'66046505640663076703'. 



Fig. 12. - Diagrams of the mean squared distances matrix, D mm i(t) for '(60)25' copolymer in 
kinetics after the same quench as in Fig. 10. Diagrams (a-d) correspond respectively to the 
following moments in time: t = 3.5, 5.75, 8.60 and 13.0. See also caption to Fig. 9 for more 
details. 



Fig. 13. - The diagram of the mean squared distances matrix, D mm i (t) for copolymer sequence 
'(60)362(60)902(60)56202(60)4' at the time moment t = 45, so that this conformation is close to 
the final equilibrium. Parameters of the quench are the same as in Fig. 10. See also caption to 
Fig. 9, 12 for more details. 
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